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1.0  INTRODUCTION 

Since  it  is  not  feasible  to  experimentally  evaluate  an  aircraft  turbine  engine  under  all 
possible  conditions,  mathematical  models  are  built  to  simulate  the  operation  of  the  major 
components  of  an  engine  such  as  the  compressor,  turbine  and  combustor  which  can  be 
tested  seperately  and  together.  In  this  report  we  are  concerned  with  the  mathematical  and 
numerical  models  of  compressors  and  turbines.  The  following  specific  topics  are  discussed: 

1.  Numerical  stability  of  a  one-dimensional  compressor  model. 

2.  Numerical  stability  of  a  multi-dimensional  finite  volume  code  that  simulates  a 
turbine. 

3.  Stage  characteristics  used  in  the  compressor  models. 

4.  Multi-dimensional  compressor  model  equations  for  the  simulation  of  circumfer¬ 
entially  distorted  flows  and  their  boundary  conditions. 

Kimzey  (Ref.  l)  developed  a  one-dimensional  time-dependent  model  for  the  analysis 
of  the  effects  of  planar  disturbances  on  a  compressor  on  the  basis  of  conservation  laws  of 
mass,  momentum  and  energy.  Experimental  data  are  used  to  synthesize  the  stage  charac¬ 
teristics  of  the  compressor.  The  compressor  stage  force  and  shaft  work  which  are  needed 
in  the  model  are  calculated  based  on  these  characteristics.  The  conservation  equations 
are  discretized  spatially  by  the  use  of  a  two-sided  difference  scheme.  Boundary  conditions 
are  imposed  on  total  pressure  and  total  temperature  at  the  inlet  boundary  and  on  static 
pressure  or  airflow  rate  at  the  exit  boundary.  The  resulting  system  of  ordinary  differential 
equations  is  integrated  forward  in  time  by  a  fourth-order  Runge-Kutta  scheme. 

These  models  have  been  applied  to  a  variety  of  compression  systems  by  Kimzey  (Ref. 
1)  and  Chamblee,  Davis  and  Kimzey  (Ref.  2).  They  are  used  to  analyze  and  extrapolate 
the  experimental  data  and  to  study  the  effects  of  unsteady  disturbances  of  the  aerodynamic 
stability  of  the  compression  system.  These  models  are  not  always  numerically  stable.  Some 
of  the  techniques  used  for  overcoming  the  numerical  instabilities  are  to  increase  the  friction 
coefficient  in  the  inlet  and  exit  ducting,  to  alter  the  ducting  lengths  or  areas  and  to  time 
average  the  numerical  solutions  once  every  few  time  steps.  However,  these  techniques  can 
achieve  only  conditional  numerical  stability  Tor  some  of  the  compression  systems.  Also  the 
application  of  the  models  has  been  restricted  to  limited  regions  of  the  operating  map.  for 
which  they  are  numerically  stable. 

Davis  (Ref.  1)  has  applied  MacCorinack's  explicit  finite  difference  scheme  to  solve  the 
partial  differential  equations  of  the  model  and  an  approximate  version  of  the  method  of 
characteristics  to  impose  the  boundary  conditions,  i  his  scheme  is  more  stable  numerically 
than  the  earlier  method,  but  this  also  exhibits  numerical  oscillations  in  the  solutions  under 
certain  conditions.  These  oscillations  are  controlled  or  reduced  by  the  addition  of  extra 
friction  or  dissipation  in  the  inlet  duct  and  by  reducing  the  inlet  duct  volume.  However. 
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the  additional  extra  friction  may  degrade  the  accuracy  of  the  simulation  of  the  actual 
physical  system.  Davis  has  applied  the  one-dimensional  compressor  model  to  simulate  the 
dynamics  of  single  spool  and  dual  spool  compression  systems.  The  computer  code  based 
on  this  scheme  is  called  COMP2SP  in  this  report  for  identification  purposes. 

Reddy  and  Tsui  (Ref.  3)  have  analyzed  the  compressor  model  equations  and  con¬ 
cluded  that  the  boundary  conditions  should  be  imposed  accurately  based  on  the  method 
of  characteristics  (M.O.C.).  They  have  tested  various  numerical  schemes  together  with 
the  M.O.C.  boundary  conditions  and  concluded  that  MacCormack  Scheme  provided  a  re¬ 
liable  and  numerically  stable  method.  They  have  not  reported  any  numerical  oscillations 
mentioned  in  Ref.  4.  The  one- dimensional  compressor  code  based  on  this  work  is  called 
COMPUT. 


One  of  the  topics  addressed  in  this  report  is  the  identification  of  the  source  of  nu¬ 
merical  oscillations  in  the  one-dimensional  two  spool  compressor  model,  COMP2SP  and 
the  changes  that  have  been  made  to  remove  the  numerical  oscillations  and  to  enhance  the 
accuracy  of  the  model.  This  is  essentially  accomplished  by  detailed  comparisons  of  the 
results  of  COMP2SP  and  COMPUT  for  an  identical  test  case,  as  reported  in  Section  2. 

The  ARO-1  code  is  a  three-dimensional,  finite  volume  cade  for  solving  the  unsteady 
Euler  equations  using  an  explicit  MacCormack  scheme.  This  code  has  been  modified  at 
AEDC  to  simulate  a  turbine  by  including  source  terms,  such  as  the  blade  force  and  shaft 
work.  The  source  terms  are  computed  from  the  known  performance  characteristics  of 
the  turbine.  This  is  called  the  ATAC  code.  Numerical  oscillations  were  observed  in  the 
solutions  of  this  code  in  and  around  the  turbine.  In  order  to  analyze  this  problem,  a  turbine 
version  of  the  COMPUT  code  has  been  developed  and  it  did  not  produce  any  numerical 
oscillations.  From  detailed  comparisons  it  has  been  found  that  the  source  of  oscillations 
in  the  ATAC  code  was  the  manner  in  which  the  source  terms  were  incorporated  into  the 
algorithm.  As  reported  in  Section  3,  appropriate  changes  have  been  made  to  the  ATAC 
code  which  remove  the  spurious  oscillations.  Some  transient  flow  results  based  on  the 
one-dimensional  turbine  code  have  also  been  reported. 

The  compressor  model  equations  typically  depend  on  various  stage  characteristics. 
These  are  obtained  experimentally  at  great  expense  and  can  be  done  only  for  some  dis¬ 
crete  rotational  speeds.  Intermediate  stage  characteristics  have  to  be  obtained  based  on 
some  interpolation  techniques.  It  is  necessary  to  develop  some  rational  models  for  vari¬ 
ous  stage  performance  characteristics  in  order  to  interpolate  for  intermediate  speeds.  The 
compressor  models  developed  at  AEDC  have  traditionally  used  the  stage  characteristics 
based  on  curves  representing  the  pressure  coefficient.  tyt,  and  temperature  coefficient,  n \ 
as  functions  of  the  flow  coefficient,  <j>  for  various  corrected  speeds.  Polynominal  surface 
fits  for  interpolating  the  stage  characteristics  for  intermediate  corrected  speeds  have  been 
found  to  be  unsatisfactory.  We  report  a  study  of  using  alternate  stage  characteristics  which 
may  be  interpolated  better  for  intermediate  speeds,  fn  Section  4.  a  calculation  procedure 
based  on  the  total  pressure  loss  coefficient,  -Si  and  the  deviation  angle,  6  as  functions  of  the 
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angle  of  incidence,  i  and  corrected  speed  is  outlined.  Based  on  test  results,  stage  charac¬ 
teristics  based  on  and  6  as  functions  of  i  are  less  sensitive  to  errors  in  interpolation  for 
intermediate  speeds  than  ifr p  and  vs  <j>. 

Multi-dimensional  compressor  models  are  necessary  for  studying  the  effects  of  dis¬ 
tortion  in  the  inlet  Bow  and  other  non-uniformities  in  the  flow  field  such  as  the  rotating 
stall  on  the  compressor  operation.  In  Section  5,  unsteady  two-dimensional  compressor 
model  (z  —  9  model)  equations  have  been  developed  and  analyzed.  This  is  a  radially 
averaged  model  suitable  for  circumferentially  distorted  flows.  Characteristics  and  compat¬ 
ibility  equations  which  can  be  used  to  impose  the  inflow  and  outflow  boundary  conditions 
accurately,  have  been  developed.  These  equations  are  developed  in  a  form  suitable  for 
developing  a  two-dimensional  compressor  code  using  similar  numerical  techniques  as  those 
used  in  one-dimensional  codes. 
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2.0  ONE-DIMENSIONAL  COMPRESSOR  MODEL 

As  discussed  in  the  last  section,  the  one  dimensional  two  spool  compressor  model, 
developed  at  AEDC,  exhibited  numerical  oscillations  under  some  conditions  and  was  un¬ 
stable  under  some  other  conditions.  For  identification  purposes  we  name  the  computer 
code  for  this  model  as  COMP2SP  in  this  report.  We  describe  in  this  section  the  changes 
that  have  been  made  to  COMP2SP  to  remove  the  numerical  oscillations  and  instability  and 
to  enhance  the  accuracy  of  the  model.  The  source  of  numerical  oscillations  in  COMP2SP 
was  identified  by  comparing  the  results  of  this  code  with  the  results  of  a  one-dimensional 
single  spool  compressor  model  developed  earlier  at  UTSI,  which  we  name  COMPUT,  for 
an  identical  test  case.  The  test  case  run  was  for  a  ten  stage  compressor  operating  at  87 
percent  corrected  speed.  A  brief  description  of  the  salient  features  of  the  two  models  is 
given  below,  followed  by  the  changes  made  to  COMP2SP  and  the  results. 

2.1  COMPUT  (SINGLE  SPOOl) 

The  governing  equations  are  derived  (Refs.  1  and  2)  by  the  application  of  the  mass, 
momentum  and  energy  conservation  laws  to  an  elemental  control  volume  in  which  the 
blade  forces,  wall  shear  forces,  shaft  work  done,  heat  added  to  the  fluid  and  mass  bleed 
How  rate  are  included.  The  resulting  system  of  first  order  partial  differential  equations  can 
be  written  as  follows: 

+  +  ^“>*’0  =°  (1) 

where, 

pA  '  pAU 

u(z,t)  =  pAU  .  f[u)  =  A(pU2  +  P) 

.pA{e- f-^)J  pAUCpTt  . 


r  wb  l 

fli 

g(u,3,t)=  -P^-F 

— 

92 

(2) 

~WS  -  Q  J 

-53- 

The  various  symbols  represent  the  following: 
p  -  density 
A  -  area 

U  -  axial  velocity 
e  -  internal  energy 
P  -  static  pressure 

Cp  -  specific  heat  at  constant  pressure 
Tt  -  stagnation  temperature 
WB  -  compressor  bleed  flow  rate 

F  -  force  of  compressor  blading  and  casing  friction  acting  on  the  fluid 
iVj  -  stage  shaft  work  added  to  fluid  in  control  volume 
Q  -  rate  of  heat  addition  to  control  volume 
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In  conjunction  with  the  partial  differential  Eq.  (1),  we  have  the  equation  of  state  given  by, 


P  —  pRT 


{3) 


where  T  is  the  static  temperature  and  It  the  gas  constant.  In  addition  we  have  the 
equations  relating  the  internal  energy  and  stagnation  temperature  to  other  variables  as 
follows: 

e  =  CpTh  (4) 

CPTt  =  CpT+ ^  (5) 

The  area  distribution  A(z)  is  known  for  a  given  compressor  system.  In  the  vector  g, 
which  acts  like  a  forcing  function  in  the  differential  equation,  various  terms  are  modeled 
empirically  for  a  particular  compressor.  F(z,  t )  represents  the  forces  of  the  blades  and  the 
casing  friction  acting  on  the  fluid.  In  practice  it  is  difficult  to  isolate  F(z,  f)  empirically 
from  the  experimental  data  and  hence  the  total  term  (F  +  P|~)  which  represents  the  forces 
including  the  wall  pressure  area  force  is  modeled  from  the  experimental  data.  Wa[z,t)  is 
the  shaft  work  done  on  the  fluid,  which  is  calculated  from  the  stage  characteristics  of 
the  compressor.  The  stage  characteristics  are  modeled  empirically  based  on  experimental 
measurements  of  stage  total  temperature,  flow  rate,  total  pressure  and  flow  angularity  at 
the  stage  entry  and  exit. 

In  Ref.  3,  various  numerical  schemes  for  integrating  these  equations  and  their  sta¬ 
bility  characteristics  have  been  analyzed.  It  has  been  concluded  that  McCormack  scheme 
together  with  the  method  of  characteristics  (M.O.C)  boundary  conditions  is  more  accurate 
and  reliable  than  other  schemes  tried  in  predicting  the  compressor  surge.  With  this  scheme 
no  numerical  instabilities  have  been  encountered. 


In  Ref.  3.  the  following  compatibility  equations  along  the  three  characteristics  of  the 
system  of  Eq.  (i)  have  been  derived. 


dP  2. 

dT+c*' 


dz 

h  =  0  along  —  =  V 


(6) 


dU  dP  _  .  ,  dz 

pc-—  ~  Pc9z  ~  9z  =  0  along  ~n  =  v  *  c 
dt  at  at 

dU  dP  ,  dz  r, 

where,  p  —  pA.  P  -  PA  and 


'y  i  ‘ 

r 

/  i 

* 

92 

- 

,T!/l  "  =,92 

■  9z- 

-  fi  -  i)u9i  -(• 

7  -  l )  fl.i 

(7) 

(8) 


(9) 
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where  gi,gz  and  g$  are  given  by  Eq.  (2).  For  subsonic  flows,  two  boundary  conditions  are 
imposed  at  the  left  boundary.  Typically,  the  total  pressure  P^(f)  and  total  temperature 
Tt(t)  are  prescribed  at  the  inflow  boundary.  These  two  boundary  conditions  and  the 
compatibility  Eq.(2.8)  are  salved  to  obtain  all  the  flow  variables  at  the  inflow  boundary, 
which  are  necessary  for  the  interior  point  numerical  scheme.  Similarly,  one  boundary 
condition,  say,  the  prescription  of  the  static  pressure  P(t)  at  the  outflow  boundary  and 
the  two  compatibility  Eqs.  (6)  and  (7)  along  the  out-going  characteristics  are  sufficient  to 
compute  all  the  variables  at  the  outflow  boundary.  Details  of  this  scheme  are  given  in  Ref. 
3. 


2.2  COMP2SP  (TWO  SPOOL) 

The  one-dimensional  compressor  modeling  technique  has  been  extended  to  include 
dual-spool  compression  systems  of  the  type  shown  in  Fig.  1.  and  has  been  applied  to  a  13- 
stage  compression  system  by  Davis  in  Ref.  4.  Variable  compressor  inlet  geometry  is  used 
on  both  the  fan  and  high  pressure  compressor  and  is  scheduled  as  a  function  of  airflow, 
which  is  determined  from  compressor  inlet  temperature  and  compressor  rotor  speed.  Stage 
characteristics  necessary  for  model  construction  were  based  on  experimental  data.  High 
pressure  compressor  characteristics  were  the  same  as  those  used  in  the  single  spool  model. 
Fan  characteristics  were  synthesized  from  the  fan  rig  data. 


COMP2SP  differs  mainly  from  COMPUT,  in  the  boundary  treatment.  In  Ref.  4, 
ignoring  the  viscous  forces  and  conduction  heat  transfer  along  the  characteristic  curves, 
the  compatibility  equations  were  derived  to  be, 


dP  -  pcdU  =  0  for 


dz 

dt 


dP  +■  pcdU  —  0  for 


dz 

dt 


=  U  -  c 

(10) 

-  U  +  c 

(U) 

Inlet  Boundary  Solution:  The  inlet  thermodynamic  properties  for  any  time  can  be  cal¬ 
culated  by  specifying  certain  boundary  conditions  ( Pt  and  Tt)  and  using  the  characteristic 
relationships  of  Eq.  (10).  The  characteristic  equation  is  first  solved  by  approximating  the 
total  derivatives  by  differences, 


Zncw  -  Zi  ~  [U  -  c)At  ( L2a) 

Figure  2  shows  the  determination  of  the  intersection  of  the  characteristic  curve  (point  zf) 
with  the  geometry  of  the  previous  time  step  (line  AB).  The  thermodynamic  relationships 
at  point  3/  can  be  determined  by  linear  interpolation  of  the  properties  between  points 
A  and  B.  Once  this  point  is  known,  the  compatibility  Eq.  (10)  can  be  approximated  by 
differences 

Pz„.  -  Pzt  ~  pc(6'r„..>t  —  U*,)  (126) 
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An  iterative  technique  is  employed  to  solve  Eq.  (12b)  along  with  the  boundary  conditions 
that  Pt  and  Tt  are  prescribed.  To  obtain  a  more  accurate  value  of  ,  an  outer  iteration 
is  applied  to  Eq.  (10),  with  density  and  acoustic  velocity  now  being  an  average  between  the 
interpolated  value  in  the  old  time  frame  and  the  previous  iterated  value  at  the  new  time 
step.  All  thermodynamic  properties  are  calculated  using  these  relationships  and  equation 
of  state. 

Exit  boundary  solution:  For  subsonic  exit  conditions,  static  pressure  is  specified  as 
the  exit  boundary  condition.  A  characteristic  scheme  similar  to  the  one  used  at  the  inlet 
boundary  can  be  employed,  fn  addition  to  using  the  characteristic  relationships  of  Eq. 
(11),  the  third  characteristic  equation  along  a  streamline  is  also  used. 

dP-c2dp  =  0  for  ~  =  U  (13) 

at 

A  procedure  similar  to  the  inlet  solution  procedure  is  used  with  each  compatibility  equation 
solved  along  its  characteristic  curve  or  streamline  curve  as  shown  in  Figure  3.  With  the 
specification  of  static  pressure  and  the  two  compatibility  equations,  iterative  procedure 
is  not  needed  in  this  case,  ft  may  be  noted  that  the  Eqs.  (2.14),  (13)  and  (IL)  are 
approximations  of  the  more  general  Eqs.  (6),  (7)  and  (8),  respectively. 

2.3  MODIFICATIONS  MADE  TO  COMP2SP 

The  following  modifications  were  made  to  the  COMP2SP  code  to  remove  numerical 
oscillations  and  instabilities. 

1.  An  arithmetic  error  in  the  subroutine  FRICFZ.  which  computes  friction  forces  in  the 
duct  volumes,  is  corrected. 

2.  The  flow  variables  in  the  FORCES  subroutine,  which  computes  the  force  terms  for 
compressor  volumes,  are  updated  after  the  predictor  step  of  the  finite  difference 
scheme,  in  addition  to  updating  them  after  the  corrector  step.  Also,  friction  forces  in 
the  last  duct  volume  are  included  in  the  scheme  which  was  not  done  in  the  original 
code. 

3.  AceuraLc  boundary  conditions  are  imposed  based  on  the  method  of  characteristics 
at  both  inflow  and  outflow  boundaries  using  Eqs.  (6),  (7)  and  (8)  as  was  done  in 
COMPUT  code. 

With  change  l.  COMP2SP  does  not  exhibit  numerical  oscillations  and  yields  stable  calcu¬ 
lations  as  shown  in  Figs.  4  and  5  for  400  time  steps. 

Figures  6  and  7  show  the  calculations  of  COMP2SP  with  changes  1  and  2.  They  are  also 
stable  and  are  almost  identical  to  Figs.  4  and  5. 

Figures  8  and  9  show  the  calculations  of  CO.V1P2SP  with  changes  1.  2,  and  3.  They  show 
that  the  How  reaches  stable  steady  asymptotic  state  in  only  300  time  steps.  These  results 
are  slightly  more  accurate  than  those  of  Figs.  4-7.  as  can  be  expected  by  the  more  accurate 
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boundary  conditions.  With  these  changes  COMP2SP  becomes  identical  to  COMPUT  in 
a  single  spool  case. 
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3.0  TURBINE  MODELING 

A  general  three-dimensional,  finite  volume  code  that  solves  unsteady  Euler  equations 
using  an  explicit  MacCormack  scheme  has  been  modified  at  AEDC  to  simulate  a  turbine 
by  the  addition  of  appropriate  forcing  functions  and  by  the  introduction  of  proper  ge¬ 
ometry  and  boundary  conditions.  The  forcing  functions  are  calculated  based  on  turbine 
stage  characteristics  which  are  input  as  data.  This  turbine  modeling  code  is  named  the 
ATAC  code.  The  testing  of  this  code  has  been  carried  out  in  one-dimensional  mode.  In 
one-dimensional  simulation  of  a  single  stage  turbine  the  code  has  exhibited  non-physical 
numerical  oscillations.  We  have  analyzed  the  source  of  these  oscillations  by  comparing  the 
results  of  the  ATAC  code  with  the  results  of  a  turbine  version  of  our  one-dimensional  finite 
difference  code  (COMPUT),  which  we  name  TURBLIT.  The  source  of  the  spurious  oscil¬ 
lations  in  ATAC  code  has  been  identified  to  be  in  the  calculation  of  the  forcing  functions. 
This  section  describes  the  details  of  the  calculation  procedures  for  the  forcing  functions  in 
both  TURBUT  and  ATAC  codes  and  the  changes  made  in  ATAC  to  remove  the  numeri¬ 
cal  oscillations.  We  also  report  some  results  obtained  by  TURBUT  under  transient  flow 
conditions  in  a  turbine. 

The  source  terms  are  computed  as  functions  of  Mach  number  from  the  known  per¬ 
formance  characteristics  of  the  turbine.  In  the  test  case,  a  single  stage  turbine  is  divided 
into  eight  equal  volumes.  The  computational  region  is  a  long  constant  area  duct  extending 
from  -5.7  to  5.9  ft.  The  turbine  inlet  is  at  0.0  ft  and  the  outlet  is  at  0.2  ft.  The  entire  duct 
is  divided  into  80  control  volumes. 

Initially  the  flow  is  assumed  to  be  uniform  from  the  inlet  of  the  duct  to  the  inlet  face 
of  the  turbine,  and  from  the  exit  face  of  the  turbine  to  the  duct  outlet.  The  flow  through 
the  turbine  is  computed  from  the  performance  characteristics.  Starting  from  these  initial 
conditions,  the  ATAC  code  is  run  at  AEDC  for  1800  time  steps,  at  a  CFL  number  or  0.8. 
The  steady  state  distribution  of  density,  mass  flux  and  Mach  Number  with  axial  distance 
are  shown  in  Figs.  10-12.  As  can  be  seen  the  oscillations  in  these  profiles  are  localized, 
just  inside  and  near  the  turbine  inlet  and  outlet,  and  hence  could  not  be  attributed  to  any 
problems  at  the  boundaries.  Further,  the  initial  inlet  Mach  Number  to  the  turbine,  which 
is  close  to  the  final  steady  state  value,  had  to  be  held  fixed  during  the  calculations.  When 
the  inlet  Mach  number  was  allowed  to  vary  with  time,  the  code  became  unstable  in  a  few 
time  steps. 

In  order  to  solve  this  problem,  a  one-dimensional  compressor  code,  COMPUT  de¬ 
scribed  in  the  previous  section  has  been  modified  to  handle  the  turbine  performance  char¬ 
acteristics.  The  modified  code  is  called  TURBUT.  This  is  a  finite  difference  code  which 
uses  an  explicit,  predictor-corrector,  MacCormack  scheme.  When  the  code  is  run  for  4.800 
time  steps,  there  is  no  evidence  of  oscillations  in  the  flow  variables,  as  can  be  seen  from 
Figs.  I.j-15.  The  cause  of  the  oscillations  in  ATAC  code  is  traced  to  the  manner  in  which 
the  source  terms  are  handled  in  that  code. 
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3.1  FORCE  FUNCTIONS  IN  THE  FINITE  DIFFERENCE  SCHEME  (TURBUT) 

For  a  specific  turbine  volume,  knowing  the  inflow  conditions  and  the  stage  character¬ 
istics,  the  outflow  conditions  can  be  calculated.  The  differences  between  the  outflow  and 
inflow  values  of  impulse  functions  give  the  forces  and  the  shaft  work.  How  the  blade  forces 
are  included  in  the  momentum  equation  is  described  below.  Suppose  the  computation  is  at 
the  node  i.  Using  the  flow  conditions  at  (i  —  l)  station  as  inflow  variables,  outflow  variables 
at  the  end  of  a  single  turbine  volume  (station  i)  can  be  computed  using  the  turbine  stage 
characteristics.  The  force  function  at  the  (i-l)  station  is  given  by 

g2{i-l)  =  IMPtn-IMPout 


where, 

IMPin  =  impulse  function  defined  as,  A(ptt 2  +  p) 
calculated  at  station  (t  -  l) 

IMPout  —  Impulse  function  calculated  from  stage  characteristics 
at  station  i 

Consider  the  momentum  equation, 

For  point  t,  the  predictor  step  becomes, 


a  IMP{i)  —  lMP(i  -  l)  .  ft(i--l) 

a^fAU)  '• - £ - T  ~£T~ 


1 

Az 


-  tAf P{i  -  L)  -  (fMPuut 


IMPxn) 


IMP,„  =  IMP[i  -  1)  and  in  the  steady  state,  IMPout  =  IMP(i)  so  that  the  time 
derivative  is  zero.  However,  in  general  IM Pout  ^  fAfP(i) 

For  the  corrector  step  we  have, 


dt 


(MO  ;.= 


fMP{i  -r  1)  -  IMP(i)  g-i  (i) 
Az  +  Az 


l 

Az 


1MP{i 


I)  !AfP{i)  -  [IMP^t  ~  fM Pm) 
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Once  again,  IMPin  -  lMP{i)  and  for  the  asymptotic  steady  state  condition,  IMP0Ut  — 
IMP(i  +  1),  so  that  the  time  derivative  vanishes  only  in  the  asymptotic  steady  state. 

Thus  proper  switching  of  the  index  for  the  source  terms  in  the  predictor  and  corrector 
steps  gives  the  correct  steady  state  values. 

3.2  FINITE  VOLUME  SCHEME  (ATAC) 

The  finite  volume  scheme  approximates  the  governing  equations 

du  - 

in  integral  form  over  a  control  volume  V . 

I-  [  udV  +  /  V  •  fdV  +  [  gdV  =  0 

df  Jv  Jv  Jv 

The  numerical  solution  algorithm  used  to  solve  the  above  system  of  equations  is  an  explicit 
predictor-corrector  scheme. 

For  predictor, 

«,V+1  =  u'v  -f-ojAt  +  5^/iT  'sk  j  -s 

For  corrector, 

-  E*'  H  ♦  (‘  -  c)  ■'  +  f  [-V  ( £  *  ■  *‘+  +  £  A  •  .* )  - » 

Consider  the  approximation  for  the  axial  momentum  equation  applied  to  a  one-dimensional 
control  volume  Vt  for  illustrating  the  detail. 


Predictor: 


\{pAU)  |,=  -  -1-  IMP(t)  -  IMP{i-  1 .)  +*a(»‘  -  l) 


where, 


Corrector: 


<?.(l  -  I)  =  IMP, a  -  l M Pout 

lMP[n  -  f\fP(i  -  l) 

IM P.„it  -  IMP  Computed  based  on 

stage  characteristics  for  volume  z 

£(p,U:)  !  -  •  ~~  lMP{i  ±  1)  -  IM P(i)  - 
Ol  i  A: 
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where, 

g2{i)  =  IMPm  ~  IM Pout 
IMPtn  =  I  MP(i) 

IMPout  -  IMP  computed  based  on  stage 
characteristics  for  volume  (i+1) 


The  above  description  of  the  Force  calculation  in  the  corrector  step  of  the  finite  volume 
scheme  is  the  modified  version  in  the  ATAC  code,  [n  the  earlier  version  of  the  corrector, 
the  force  calculation  for  V,  was  based  on  flow  conditions  in  V,  as  inflow  conditions  for 
which  was  making  the  scheme  to  produce  spurious  numerical  oscillations. 

Subsequently,  the  corrector  setp  has  been  modified  in  such  a  way  that  the  force  on  Vt  is 
based  on  flow  conditions  at  Vt  going  into  Vt+  (,  which  makes  the  scheme  stable  and  remove 
the  spurious  oscillations.  It  may  be  noted  that  in  the  modified  version,  time  derivatives 
go  to  zero  in  the  asymptotic  steady  state  as  they  should.  Modified  ATAC  code  in  one- 
dimensional  case  now  produces  the  same  results  as  the  one-dimensional  code  TURBUT. 

3.3  STUDIES  OF  TRANSIENT  FLOW  CONDITIONS  IN  A  TURBINE 

In  order  to  study  the  effects  of  transient  flow  conditions  on  a  turbine,  the  code  TUR¬ 
BUT  was  run  with  changes  in  the  static  pressure  at  the  exit. 

Computation  was  started  with  steady  state  flow  conditions  at  the  time  step  6600, 
established  in  an  earlier  run.  At  that  time  the  static  pressure  at  the  outflow  plane  of  the 
exit  duct  was  impulsively  dropped  from  9.485  psi  to  8  psi.  By  12,600  time  steps,  i.e.,  in  an 
additional  6,000  time  steps  the  flow  reached  a  new  steady  state  condition.  The  time  step 
used  was  0.000015  sec,  with  the  Courant  number  0.8.  Figures  16  and  17  show  the  plots 
of  static  presure,  and  mass  flow,  at  various  time  steps  in  the  duct  which  has  the  turbine 
volumes  located  between  0  and  0.2  units.  It  was  noticed  that  for  several  hundred  time  steps 
most  of  the  transient  phenomena  was  limited  to  the  exit  duct,  which  eventually  propagated 
into  the  inlet  duct  as  well.  The  Mach  number  and  mass  flow  in  the  exit  duct  underwent 
strong  transient  perturbations,  and  not  as  much  in  the  inlet  duct,  before  settling  down  to 
new  steady  state  conditions  compatible  with  the  pressure  at  the  exit  plane  of  the  duct. 
The  turbine  seems  to  act  tike  a  big  damper  in  not  letting  the  large  perturbations  of  the 
exit  duct  go  upstream  into  the  inlet  duct.  It  was  also  noted  that  the  head  of  rarefaction 
wave  propagated  at  a  speed  close  to  the  propagation  speed  of  a  small  disturbance. 

The  above  experiment  was  repeated  with  a  sudden  increase  in  exit  pressure  from  9.485 
psi  to  10  psi.  As  before,  computation  was  started  with  steady  state  How  conditions,  at 
the  time  step  6600.  At  that  time  the  static  pressure  at  the  outflow  plane  of  the  exit  duct 
was  impulsively  increased  to  10  psi  and  held  fixed.  The  relative  magnitude  of  the  pressure 
pulse  was  0.054.  In  an  additional  4,500  time  steps  the  flow  had  reached  a  new  steady  state 
condition.  Figure  L8  is  a  plot  of  the  non-dimensional  static  pressure  in  the  duct  at  various 
time  steps.  The  CFL  number  for  this  computation  was  0.8. 

The  next  set  of  experiments  consisted  of  introducing  a  compression  wave  Followed  by  a 
rarefaction  wave.  This  was  accomplished  by  impulsively  increasing  the  exit  pressure  from 
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9.485  psi  to  10  psi,  holding  it  fixed  for  LOO  time  steps,  and  then  dropping  it  back  to  9.485 
psi.  Figure  19  is  a  plot  of  static  pressure  vs  distance  downstream  of  the  turbine  exit  at 
various  time  steps.  At  the  time  step  6600,  the  exit  pressure  was  increased  to  10  psi  and 
held  fixed  for  100  time  steps.  At  the  time  step  6700,  the  exit  pressure  was  dropped  back 
to  9.485  psi  and  held  fixed  thereafter.  This  set  up  a  pressure  pulse  travelling  up  the  duct. 
The  amplitude  of  the  pressure  pulse  was  0.415  psi.  Figure  20  is  similar  to  Fig.  19  with 
plots  drawn  at  closer  time  intervals.  The  wave  speed,  obtained  from  Fig.  20  by  calculating 
the  time  it  takes  for  the  peak  of  the  pulse  to  travel  from  one  location  to  another,  is  828 
fps.  This  value  agrees  well  with  the  local  ( U  -  c)  value  (acoustic  propagation  speed)  of 
838  fps  that  exists  at  the  peak  of  the  wave. 

The  leading  edge  of  the  pressure  pulse  hits  the  exit  face  of  the  turbine  at  the  time 
step  7000  (t  =  0.006  sec),  and  is  almost  instantaneously  transmitted  through  the  turbine. 
It  may  be  noted  that  the  turbine  extends  a  distance  of  only  0.2  units.  The  propagation 
of  the  wave  that  is  transmitted  upstream  through  the  turbine  is  shown  in  Fig.  21.  The 
amplitude  of  the  wave  is  0.121  psi.  While  the  amplitude  of  the  transmitted  wave  is  reduced, 
the  propagation  speed  increases  gradually  from  the  828  fps  to  the  local  acoustic  speed  of 
924  fps  in  the  upstream  duct. 

Figure  22  is  a  plot  of  the  static  pressure  versus  distance  in  the  downstream  portion 
of  the  duct  for  N  >  7000.  where  N  is  the  number  of  time  steps.  The  effect  of  the  multiple 
reflections  of  the  pulse  at  the  turbine  exit  face  and  the  exit  of  the  duct  on  the  pressure 
distribution  can  be  clearly  seen.  As  can  be  seen  from  the  N  =  7400  values,  the  flow  will 
settle  to  the  correct  steady  state  value. 
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4.0  STAGE  CHARACTERISTICS 

The  compressor  model  equations  typically  depend  on  various  stage  characteristics.  These 
stage  performance  characteristics  are  obtained  experimentally  at  great  expense  and  can  be 
done  only  for  some  discrete  rotational  speeds.  Intermediate  stage  characteristics  have  to 
be  obtained  based  on  some  interpolation  techniques.  Simple  polynomial  surface  fits  have 
been  found  to  be  unsatisfactory  for  stage  characteristics  based  on  pressure  coefficient  and 
temperature  coefficient  versus  flow  coefficient.  It  is  necessary  to  investigate  alternate  stage 
performance  characteristics  suitable  for  high  speed  compressors  in  order  to  interpolate  for 
intermediate  speeds. 

4.1  CALCULATION  OF  STAGE  CHARACTERISTICS 


In  lieu  of  the  parameter  flow  coefficient  <j> ,  which  is  normally  used,  the  calculation 
procedure  proposed  mainly  uses  two  parameters,  namely,  total  pressure  loss  coefficient  (w) 
and  deviation  angle  (5),  represented  as  functions  of  two  variables  -  corrected  speed  and 
angle  of  incidence.  The  deviation  angle  and  angle  of  incidence  are  defined  as  shown  in 
Fig.  23.  In  addition  to  their  common  use,  the  choice  of  w  and  6  is  preferable  to  other 
parameters,  since  both  are  relatively  insensitive  to  Mach  numbers  or  corrected  speeds  when 
corrected  speed  is  low  and  also  not  too  sensitive  to  angle  of  incidence.  As  a  test  case  the 
eight-stage  axial  flow  compressor  of  the  J85-I3  engine  has  been  chosen  for  this  purpose. 

In  addition  to  w  and  6.  the  procedure  uses  the  stage  temperature  rise  as  calculated 
from  the  Euler  turbine  equation. 


(>  + 


*§" 


where, 

M  : 

C 

c:  acoustic  speed  at  station  in  question  e.g.  Af ,  -  ^ 

Mr  :  — 

<*'1 

Ma  :  f 

w:  axial  velocity 

N:  rpm  of  engine,  fl  =  2^ 

Tt :  total  temperature 
T:  static  temperature 
Pt-  total  pressure 

P'i-  total  pressure  with  respect  to  rotor  coordinates 
.4.  =  annulus  area 
p:  static  pressure 
p  —  density 
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Subscript  0  corresponds  to  the  variables  in  front  of  the  stator,  subscript  1  corresponds 
to  the  exit  of  stator  or  in  front  of  rotor  and  subscript  2  is  For  the  exit  station  of  the  rotor. 


*  (/^2  blade  + 


Superscript  ’  refers  to  variables  with  respect  to  rotor  coordinates. 

The  total  pressure  loss  coefficients  for  stator  and  rotor  are  defined  as  follows. 


For  stator. 


For  rotor, 


Wr 


K  1  -  fg 

PI  1  -  Pi 


Once  the  total  pressure  loss  coefficients  are  available,  the  pressure  rise  can  be  computed 
in  the  following  way.  For  stator, 


w.  = 


(l+^Af02)^  -1 


(15) 


from  which  Pi  /P0  can  be  calculated. 


For  rotor, 


M[2)^  -  %(1  + 


from  which  Pi/ Pi  can  be  calculated. 


(16) 


For  the  rotor,  using  the  geometrical  parameters.  Eqs.  (L5)  and  (16)  are  iteratively 
solved  along  with  mass  balance  equation,  piiOiAi  =  piA 1W2,  in  order  to  obtain  Tn,Ti,  ^ 
and  pi.  This  would  determine  the  rotor  exit  velocity  triangle  completely.  In  the  above 
iteration  procedure  Ti  is  obtained  from. 


Ta  =  'Mi  +  ^m2-) 

Pi  1  T 1  ( t  — 

For  the  stator,  using  the  geometrical  parameters,  fc’q.  (15)  and  mass  balance  equation 
pi\Wr,Ao  —  p\W\A\  are  iteratively  solved  to  obtain  ^-Pi  and  T|.  This  would  help  define 
the  stator  exit  velocity  triangle  completely.  In  the  above  iteration  procedure  Ti  is  obtained 
from, 

r,(i =  Mi  1  ~--W?) 
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4.2  INTERPOLATION  ROUTINE 

The  basic  input  curves  were  available  for  four  corrected  speeds,  namely,  80  percent, 
87  percent,  94  percent,  and  100  percent.  Stage  characteristics  were  determined  at  all  the 
corrected  speeds  using  the  above  mentioned  calculation  procedure.  In  order  to  check  the 
accuracy  of  the  procedure,  it  was  decided  to  generate  the  stage  characteristics  at  87  percent 
corrected  speed  by  using  the  values  at  80  and  94  percent  corrected  speeds,  by  means  of 
a  suitable  interpolation  routine.  However,  no  noticeable  pattern  in  the  basic  input  curves 
was  observed  and  hence  linear  interpolation  was  resorted  to,  between  80  and  94  percent 
corrected  speed  curves  in  order  to  generate  the  corresponding  curves  at  87  percent  corrected 
speed.  The  curves  thus  obtained  have  been  compared  to  the  actual  curves  in  Figs.  24  -  27. 

4.3  RESULTS 

Figures  24  and  25  show  the  variation  of  incidence  angle  with  deviation  angle  for  rotors  1 
and  4,  and  correspondingly  Figs.  26  and  27  show  the  variation  of  incidence  angle  with  total 
pressure  loss  coefficient  for  both  actual  and  interpolated  values.  In  order  to  ascertain  the 
effect  of  linear  interpolation,  the  stage  characteristics  obtained  using  the  new  calculation 
procedure  for  both  actual  and  interpolated  cases,  have  been  fed  to  an  available  code  of 
a  one-dimensional  time-dependent  compressor  model,  COMPUT  and  some  representative 
values  have  been  compared  as  the  code  reaches  steady  state. 

Figure  28  shows  the  variation  of  inlet  mass  flow  with  time  as  the  steady  state  is 
reached  for  both  actual  and  interpolated  cases  and  similarly  the  variation  of  outflow  total 
temperature  with  time  is  shown  in  Fig.  29.  As  can  be  seen  from  Fig.  28  and  29,  the  inlet 
mass  flow  differs  from  the  actual  value  by  0.61  percent  and  the  total  temperature  rise  across 
the  compressor  by  l.65percenf.  The  last  value  is  close  to  the  largest  deviation  observed  in 
this  run. 

In  order  to  compare  the  method  against  conventional  way  of  feeding  the  stage  char¬ 
acteristics  to  the  code  (essentially,  How  coefficient  <j>  vs  pressure  coefficient  i L'p  and  flow 
coefficient  <j>  vs  temperature  rise  coefficient  ip?},  the  code  was  fed  with  stage  characteristics 
of  6  vs  i bp  and  (j>  vs  and  run  until  the  steady  state  is  reached.  Figure  30  shows  the 
variation  of  inlet  mass  flow  with  time  as  the  steady  state  is  reached,  and  similarly  Fig.  3L 
shows  the  variation  of  outflow  total  temperature  with  time.  As  can  be  seen  from  the  results 
for  this  method,  while  the  variation  in  total  temperature  rise  remains  the  same  as  before, 
the  variation  in  inlet  mass  flow  is  as  much  as  3.59  percent  compared  to  the  actual  value. 


*lhus  it  has  been  noticed  that  in  the  worst  case  where  there  is  no  noticeable  pattern 
in  the  input  curves  and  linear  interpolation  is  resorted  to.  the  results  are  within  2  percent 
of  the  actual  values  by  the  method  of  stage  characteristics  outlined  in  this  section.  The 
error  margin  is  expected  to  come  down  if  the  input  curves  are  more  amenable  to  a  betler 
interpolation  routine.  In  any  case,  the  current  method  is  more  suitable  for  high-speed 
compressors  than  using  \pf,  and  ky  vs  o. 
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5.0  MULTI-DIMENSIONAL  COMPRESSOR  MODELING 

For  studying  the  effects  of  distortion  in  the  inlet  flow  and  other  non-uniformities  in 
the  flow  field  such  as  rotating  stall  on  the  compressor  operation,  one  needs  to  consider 
multi-dimensional  models.  Kimzey  and  his  associates  (Refs,  t  and  2)  generalized  the 
one-dimensional  model  to  account  for  certain  three-dimensional  features  and  yet  keep  the 
partial  differential  equations  one-dimensional  in  nature.  Equations  representing  the  laws 
of  conservation  of  mass,  axial  momentum  and  energy  are  written  for  each  control  volume. 
Radial  and  circumferential  momentum  fluxes  are  ignored.  Force  and  shaft  work  terms 
are  determined  from  the  empirical  steady-state  characteristics  modified  for  the  unsteady 
aerodynamic  response  of  the  blades  based  on  Goethert-Reddy  analysis  (Ref.  tl).  The 
radial  work  distribution  and  circumferential  crossflow  contribution  to  the  swirl  velocity  are 
also  accounted  for,  by  empirical  modifications  to  the  stage  characteristics.  The  differential 
equations  are  similar  to  those  used  in  parallel  compressor  models.  These  equations  are 
integrated  in  time  by  a  Runge-Kutta  method.  The  resulting  code,  while  operational  under 
some  conditions,  suffers  From  numerical  instabilities  under  other  conditions  partly  due  to 
the  numerical  integration  scheme  and  partly  due  to  the  method  of  imposing  the  boundary 
conditions. 

Tesch  and  Steenken  (Refs.  7  and  8)  and  Hosney  and  Steenken  (Ref.  10)  at  General 
Electric  have  also  considered  generalizations  of  one-dimensional  compressor  models  to  in¬ 
clude  multi-dimensional  effects  for  various  compression  systems.  Tesch  and  Steenken  have 
used  the  parallel-compressor  analysis,  which  simulates  the  compressor  by  segmenting  it 
circumferentially  into  several  one-ditnensionai  parallel  compressor  models.  Their  analysis 
includes  circumferential  total  pressure  and  total  temperature  distortions  and  circumferen¬ 
tial  redistribution  of  distorted  flows  in  blade-free  volumes.  The  equations  are  solved  by  an 
explicit  time-marching  technique.  Hosney  and  Steenken  have  developed  a  two-dimensional 
compressor  mode!  (z-r  model)  which  allows  for  radial  redistribution  of  the  flow.  This  model 
is  developed  to  account  for  dynamic  interactions  between  the  fan  and  the  compressor  in 
a  turbofan  compression  system.  For  such  a  multi-spool  system,  circumferential  variation 
of  the  flow  is  neglected  and  conservation  equations  for  mass,  axial  and  radial  momentums 
and  energy  are  written  for  a  control  volume  in  discrete  form.  Various  forcing  functions 
are  modeled  empirically  through  stage  characteristics.  Equations  are  integrated  in  time 
explicitly.  The  GE  models  have  been  applied  to  study  the  stability  and  frequency  response 
of  various  compression  systems. 

Circumferentially  distorted  flow  fields  need  to  be  studied  under  different  operating 
conditions  of  compressors.  Compressor  surge  and  stall  are  intimately  connected  with 
distorted  flow  fields.  In  this  section  we  develop  the  partial  differential  equations  for  a 
two-dimensional  compressor  model  {/.-ti  model)  In  which  flow  variations  in  axial  and  cir¬ 
cumferential  directions  are  considered  but  How  variations  in  the  radial  direction  are  aver¬ 
aged  out.  The  equations  are  analyzed  and  are  transformed  into  compatibility  equations  in 
characteristics  directions,  which  provide  an  accurate  means  of  imposing  inflow  and  outflow 
boundary  ronditions. 
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5.1  TWO-DIMENSIONAL  (z-0)  COMPRESSOR  MODEL  EQUATIONS 

We  start  with  the  Euler  equations  of  motion  for  inviscid  flow  in  conservation  form 
in  cylindrical  polar  coordinates.  These  equations  are  integrated  in  the  radial  direction 
between  the  hub  r\  and  tip  r2.  Let  u,  v  and  w  be  the  components  of  velocity  in  the  radial, 
circumferential  and  axial  directions  respectively.  It  has  been  assumed  that  the  radial 
velocity  is  much  smaller  than  the  circumferential  and  axial  velocities,  i.e.,  u  <<  v,w. 
Defining  suitable  averages  of  the  flow  field  variables,  the  radially  averaged  equations  of 
motion  can  be  written  as  follows: 

Continuity 


dp 

dt 


0  3 

fLr  +  ~[fpwLr)  +  qq(pvLt)  +  9\ 


=  0 


(17) 


where 


rpw 


pv 


Cr  ~  f  dr  =  r2  -  r, 
fT 3  rdr  i 

J  *  I 

f  prdr 

P  ~  jVdT 
J  prwdr 


f  dr 
/  pvdr 
f  dr 


ffi  = 


dr 


dr 


rpu  -  rpw~  -  pv  — 
dz  dO 


It  may  be  noted  that  while  pr 


—  pr  ,  in  general  fg  f(j 


Axial  Momentum 


3  /  -  -i-i  d 


(pirwr  -  pr)Lr 


HO 


{pvwLT)  -rg2  -  0 


US) 
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where 


w 


f  prwdr 
f  prdr 


prw 

pr 


pwwr  = 


pvw  = 


f  pw2rdr 
f  dr 
j  pvwdr 
f  dr 


puwr  -  pvw 


dr 

30 


-  ( pw 2  +  p)r 


dr 

dz 


Tn 


^1 


F  zLr 


Circumferential  Momentum: 


Fz  = 


jFzdr 

J  dr 


St 


(pv}rLr 


3  _ ,  ,  3 

—  lfiun:rLr)  +  - 


[pw  +  p)Lr 


+  g-i--0 


(19) 


where 


/  prvdr 
J  prdr 


prv 

pr 


J pv2dr 
f  dr 


pwvr  — 


j  pwvr  dr 
J  dr 


y-j  = 


puvr  - 


3r 

36 


3r 

pwvr  --- 
dz 


F„  Lr 


F* 


f  Fjdr 
j  dr 
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Energy: 

^  ^  A 

0-t(pe)rLr  +  —  {pHrw  Lr )  +  —  Lr)  +  g4  =  0 

where 

total  enthalpy,  H  =  e  +  £ 


_  _  per  _  J"  perdr 
pr  prdr 


^  =  IpHrmb 

[dr 


pH  v  = 


/  pHvdr 

~JdT 


04  = 


pH ru  -  pHv—  -  pHrw 

do  dz 


QLr-W4LT 


(20) 


Q  = 


J  Qdr 

J  dr 


=  heat  added 


777-  f  W*dr 

W.i  —  — jr— —  =  shaft  work 
J  dr 

There  are  a  equations  (continuity,  2  momentum,  energy  and  equation  of  state)  and  15  un¬ 
knowns  (ftptw1vti,p,rlw,^i,pwir.^mstfiS^r,pv3_;pHTw,^Hv).  f  is  a  geometrically 
known  quantity.  Now  we  make  an  assumption  that  fg  =  fg.  This  closure  assumption 
gives  9  additional  equations.  Then  we  have  5  unknowns  and  5  equations. 

Among  the  forcing  function  terms  gu..gA  in  the  above  equations.  g:  =  0  if  there  is 
no  mass  bleed  and  02,03,04  have  to  be  determined  empirically  based  on  suitable  data 
correlations  for  stage  characteristics  for  a  particular  compression  system. 


The  above  equations  are  of  the  form, 


,  df{  dfi 

lk~r!h  +  M+0  =  ii 


where. 


/  r Lrp  \ 

(  31  \ 

a  - 

rLrpw 

fLrpr 

•  U  - 

03 

U‘> 

fLrpK  t 

'  04 

[‘i  I) 
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/  rpw  Lr  \ 

/  pv  Lr  \ 

( pwwr  +  pr )  Lr 

fl  — 

pvw  Lr 

pwvr  Lr 

J  2  — 

[pv V  +  P )  Lr 

V  pHrw  Lr  / 

V  pHv  Lr  / 

These  equations  are  hyperbolic  which  means  that  if  we  deiine  the  Jacobian  matrices, 


Mi 


du  ’ 


dh 

du 


(22) 


The  matrix  P  =  fciM,  -  *2^2*  where  fcj  and  fc2  are  arbitrary  parameters  such  that 
fc?  +  fc?  =  1.  has  real  eigenvalues  and  a  complete  set  of  eigenvectors.  The  eigenvalues  are 

1  A 


q,q,  q  +  c  and  q  -  c 


where  q  =  klv  +  k2w,  and  c  is  the  speed  of  sound.  For  axial  compressor  flows  where  the  flow 
is  primarily  in  the  axial  direction,  the  eigenvalues  of  the  matrix  Mi  determine  the  boundary 
condition  requirements  and  procedures.  Matrix  A  has  eigenvalues:  »,»,  w  +  c  and  w  -  c. 
For  subsonic  inflows  this  would  indicate  that  we  must  prescribe  three  boundary  conditions 
(equal  to  the  number  of  positive  eigenvalues)  at  the  inflow  boundary.  We  can  prescribe  the 
total  pressure,  total  temperature  and  circumferential  velocity  or  flow  angle  at  the  inlet.  At 
the  outflow  boundary,  if  the  flow  is  subsonic,  there  is  one  negative  eigenvalue  and  we  need 
to  prescribe  one  boundary  condition.  This  is  typically  the  static  pressure  at  the  outflow. 

Equation  (21)  can  be  integrated  by  a  numerical  scheme  such  as  the  MacCormack 
scheme  as  was  done  in  the  one-dimensional  model.  The  boundary  conditions  should  be 
imposed  at  the  boundaries  by  using  the  characteristic  compatibility  equations,  which  can 
be  derived  by  suitable  transformation  of  Eqs.  (21).  These  are  given  by 


dp 

dt 


dw 

dt 


1  dp  v  ,  dp  _  J  *?£  \  - 

Jdt  +  r'd9  c2d8]  r2  38  '  9'  e2 

along  —  w  for  a  fixed  0 
dt 

dv  I  ,  dv  1  c)p.  c2  dr 
dt  rd8  -°‘ 


along  -  w  for  a  fixed  0 
dt 

l  dp  v  dw  c dc  v  dp  l  . 

- - -  4.  -  - - J-  — - 4-  <72  —  — <74  •—  U, 

pc  dt  r  d8  r  d$  per  68  pc 

along  —  =  h;  >-  c  for  a  fixed  8 
dt. 


(23) 


(24) 


(25) 
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dtv  1  dp  v  dw  c  dv  v  dp  _  t 
~dt  -  Jcdi  +  ri  “  rde~  fc?de  +®2  “  =  °’ 

dz 

along  —  =  w  -  c  For  a  fixed  6  (26) 

Details  of  the  derivation  of  the  eigenvalues  of  My  and  the  characteristic  compatibility  equa¬ 
tions  are  outlined  in  Appendix  A.  Discretization  of  the  equations  described  in  this  section 
and  the  corresponding  Computer  Code  can  be  developed  similar  to  the  one-dimensional 
compressor  model  code  COMPUT. 
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6.0  CONCLUSIONS 

A.  The  numerical  oscillations  encountered  in  a  one-dimensional  compressor  model 
(COMP2SP)  have  been  traced  to  some  parts  of  the  algorithm.  With  appropriate 
changes,  the  code  is  now  numerically  stable  and  also  it  has  been  made  more  accurate 
by  imposing  the  boundary  conditions  by  the  method  of  characteristics. 

B.  Numerical  oscillations  encountered  in  the  turbine  model  (ATAC)  are  due  to  the  man¬ 
ner  in  which  the  source  functions  are  incorporated  in  the  finite  volume  algorithm. 
With  appropriate  changes  in  the  algorithm  spurious  oscillations  have  been  eliminated 
in  this  model. 

C.  Stage  characteristics  based  on  the  total  pressure  loss  coefficient,  Q  and  the  deviation 
angle,  6  as  functions  of  the  angle  of  incidence,  i  and  corrected  speed  have  been  found 
to  be  less  sensitive  to  errors  in  interpolation  for  intermediate  speeds  than  pressure 
coefficient,  and  temperature  coefficient,  i pt  as  functions  of  the  flow  coefficient 

D.  Radially  averaged,  unsteady  two-dimensional  compressor  model  (z-0  model)  equations 
suitable  for  circumferentially  distorted  flows  have  been  derived.  Characteristics  and 
compatibility  equations  are  derived  for  imposing  the  inflow  and  outflow  boundary 
conditions. 
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Figure  1.  Schematic  of  the  13-stage,  dual-spool  compression  system  model. 
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Figure  2.  Schematic  of  inlet  characteristic  boundary  scheme. 
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Figure  3.  Schematic  of  exit  characteristic  boundary  scheme. 
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Figure  4.  Inlet  mass  flow  versus  time  with  change  1  (COMP2SP). 
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Figure  5.  Outflow  total  temperature  versus  time  with  change  1  (COMP2SP). 
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Figure  6.  Inlet  mass  flow  versus  lime  with  changes  1  and  2  (COMP2SP). 
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Figure  7.  Outlfow  total  temperature  versus  time  with  changes  1  and  2  (COMP2SP). 
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Figure  1 1 .  Steady-state  distribution  of  mass  flux  with  oscillations  (4,800  time  steps). 
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Figure  12.  Steady-state  distribution  of  Mach  number  with  oscillations  (4,800  time  steps). 
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Figure  13.  Steady-state  distribution  of  density  (TURBUT). 
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Figure  15.  Steady -state  distribution  of  Mach  number  (TURBIJT). 
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Figure  16.  Evolution  of  mass  flow  from  une  steady-state  condition  to  another  as  exit  pressure 
is  dropped  impulsively  from  9.485  psi  to  8  psi  att  =  0.0. 
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Figure  17.  Evolution  of  static  pressure  from  one  steady-state  condition  to  another  as  exit 
pressure  is  dropped  impulsively  from  9.485  psi  to  8  psi  at  t  =  0.0. 
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Figure  1 H  Evolution  of  static  pressure  from  one  steady-state  condition  to  another  us  the 
exit  pressure  is  increased  impulsively  from  9.485  to  10  psi  at  t  =  0.0. 
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H  igure  21.  I  Jpstream  prupagatiun  ol  a  pressure  pulse  in  the  upstream  side  duct  of  a  turbine. 
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Figure  22.  Interaction  of  the  pressure  pulse  with  reflections  from  the  turbine  and 
the  downstream  duct  exit. 
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Figure  4.1  Definition  of  angle  of  incidence 
and  deviation  angle 
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Figure  23.  Definition  of  angle  of  incidence  and  deviation  angle. 
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Figure  24.  Actual  and  interpolated  curves  of  deviation  angle  versus  angle 
of  incidence  (Rotor  1). 
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Figure  25.  Actual  and  interpolated  curves  of  deviation  angle  versus  angle 
of  incidence  (Rotor  4). 
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Figure  26.  Actual  and  interpolated  curves  of  total  pressure  loss  coefficient 
versus  angle  of  incidence  (Rotor  1 ). 
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Figure  27.  Actual  and  interpolated  curves  of  total  pressure  loss  coefficient 
versus  angle  of  incidence  (Rotor  4). 
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Figure  29,  Aetna]  and  interpolated  (i  vs  8  and  i  vs  to)  curves  of  outflow  total 
temperature  with  time  as  the  steady  state  is  reached. 
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Figure  30.  Actual  and  interpolated  (8  and  w  vs  i)  and  (ipp  and  ipi  vs  $)  curves  of  inlet 
mass  flow  with  time  as  the  steady  state  is  reached. 
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Figure  31 .  Aclual  and  interpolated  (ti  and  w  vs  i)  and  (ipp  and  ipt  vs  $)  curves  of  outflow  total 
temperature  with  time  as  the  steady  state  is  reached 
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APPENDIX  A 

CHARACTERISTICS  AND  COMPATIBILITY  EQUATIONS 
Equation  (21)  can  be  rewritten  as, 

du  du  dfi  du  d/2  du  ^  d/2  dr  _  ^ 

du  dt  ^  du  dz  du  36  ^  dr  36  ^ 


Substituting  for  fjj,  and  ^  we  obtain, 

t( 


and. 


dp  ,  dp  dw  v  ,  v  dp  pdv,  pv  dr 


dp  _dws  f  2  dp 


t  vr  ,  . 

\w - r  P 

X  dt  dt 


)-( 


r  d&  r  d9J  r 2  dO 

dw  dp ,  ,  vw  dp  pv  dw 

W  dz  ^2pW~dz  ~  ds'  +  '~r"dO  *  ~r  d0 

dw  du,  _vw  dr 

H - ^7)  ~  P~ 2~^Z  +Si2  —  0 

rdv  rz  do 


,  do  dv  ,  ,  dp  dw  dv,  /V2  dp  2pv  dv  L  dp, 

<■ vd~t  ~  '*)  ~  wvTz  +  p'vd7  ’fiwdJ  +  {7¥o  +  ~To  + 


dz 
,pv 2 


,pu  p  ,  dr 

"(72”  +  7^d0  +"~ 


( 

dw 

dz 


\T_dp 
2  dt 


+  pw 
dv 


+  pwv 


i)  dv 

l 

dp s 

-  V 2  dp 

t  dt 

7  -  1 

dr 

*•  2  dz 

7 w  dp- 

rLY2 

dp 

pvw  dw 

7  -  L  dz : 

r  .  ^ 

L  2r 

do 

r  d8 
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IP 


'IP 


dv  7  v  dp-  ,  7  vp  /)»„ j,* 


_L_  _  _i:'  _/ _ _ _ lH  ^  —  ly 2) 

d0  7  -  L  r  38  S  -  i  r2  r2  2  d0 

The  above  equations  can  be  written  as, 

du  du  du  dr 


^7  +  ^  -  0 


or. 


where, 


dti 

<?7 


dz 


do 


00 


L  = 


du 

d~d 


(  l  0  0  0  \ 

w  p  0  0 

!■  0  p  0 

\  vr  piv  pi;  J 
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Mi  = 


to 

P 

u 

P 


/  w 
w2 
wv 

.  wV - 
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For  the  Mi  matrix,  the  characteristic  equation  is  given  by  j  Mi  -  A/  1=  0 


|  Mi  -  A  /  i  = 


j  w  -  A  p  0 
|  0  w  —  A  0 

1  0  0  it)  -  A 

0  7  p  0  w  —  A 


0 

I 

p 

0 


=  0 


Solving  the  above  we  obtain  the  roots  as. 

A |  —  tn.  Aj  —  w,  A 3  =  iv  —  c,  A 4  —•  w  —  c 


where  ~ipj  p  -  cl 

The  matrix  T  formed  by  the  above  eigen  vectors  is. 

/I  0  l  lx 

~  0  0  4  =JS 

1  —  i'  *■ 
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’0  0  C-  c2 
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multiplying  by  T-1, 
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substituting  the  above,  the  compatibility  equations  are  derived  as. 
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or 


or 


or 


or 


dp 

dt 
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NOMENCLATURE 

A  Area 

p  density 

t  internal  energy 

P  static  pressure,  matrix  as  defined  in  Chapter  5 

Pt  total  pressure 

cp  specific  heat  at  constant  pressure 

cv  specific  heat  at  constant  volume 

•7  ratio  of  specific  heats 

T  static  temperature,  matrix  as  defined  in  Appendix  A 
Tt  total  temperature 

Wg  compressor  bleed  flow  rate 

t  time 

F  force  of  compressor  blading  and  casing  friction  acting  on  fluid 
Wa  stage  shaft  work  added  to  fluid  in  control  volume 

R  gas  constant 

x ,y,z  coordinates  in  the  cartesian  coordinate 

system 

c  speed  of  sound 

CN  courant  number 


a  as  defined  in  Eq.  (2) 

/  as  defined  in  Eq.  (2) 

g  as  defined  in  Eq.  (2),  and  Eq.  (21) 

IMP  impulse  function 
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N  number  of  time  steps,  rpm  of  engine 

<j>  flow  coefficient 

i/ip  pressure  coefficient 

\frt  temperature  coefficient 

i  angle  of  incidence 

6  deviation  angle 

V  absolute  velocity,  volume,  total  velocity 


V'  velocity 

r  radius  at  pitchline 

1 3  absolute  flow  angle 

/?'  relative  flow  angle 

Q  total  pressure  loss  coefficient 

Ma  axial  Mach  number 

Mt  rotor  Mach  number  at  pitch  radius 

rt  hub  radius 

r-2  tip  radius 

r.  6,  z  cylindrical  polar  coordinates 

u,v,  tu  components  of  velocity  in  the  radial,  circumferential,  and  axial  directions,  respectively 
fi,fi  as  defined  in  Eq.  (21) 

A/,,  M2  Jacobian  matrices  as  defined  in  Eq.  (22),  matrices  as  defined  in  Appendix  A 
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n 


L 

M 

H 

NUN^NZ 

A|,2sa,4 

9 

Q 

SUBSCRIPTS 

V 

‘s’ 

‘O' 


V 


‘2’ 


'8' 


V 

SUPERSCRIPTS 


rotor  angular  velocity 
averaged  values  in  radial  direction 
length,  matrix  as  defined  in  Appendix  A 
Mach  number 
total  enthalpy 

matrices  as  defined  in  Appendix  A 
roots  of  the  characteristic  equation 
as  defined  in  Section  5 
rate  of  heat  addition 
arbitrary  constants 

refers  to  rotor,  radial  direction 

refers  to  stator 

in  front  of  the  stator 

after  the  stator  and  in  front  of  the  rotor 

after  the  rotor 

refers  to  circumferential  direction 
refers  to  axial  direction 


’  variables  with  respect  to  rotor  coordinates 

+-  refers  to  values  at  the  forward  node 

values  at  the  previous  node 
~  as  defined  in  Appendix  A 


68 


